On the robustness of scale invariance in SOC models 
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A random-neighbor extremal stick-slip model is intro- 
duced. In the thermodynamic limit, the distribution of states 
has a simple analytical form and the mean avalanche size, 
as a function of the coupling parameter, is exactly calculable. 
The system is critical only at a special point J c in coupling pa- 
rameter space. However, the critical region around this point, 
where approximate scale invariance holds, is very large, sug- 
gesting a mechanism for explaining the ubiquity of power laws 
in Nature. 

PACS number(s): 05.40. +j, 05.70.Ln, 64.60.Lx, 91.30.Bi. 



I. INTRODUCTION 

Self-organized criticality (SOC) is an intriguing con- 
cept which started a large 'avalanche' of research on 
mechanisms leading to scale invariance in extended dy- 
namical systems |IJ. However, there is no general agree- 
ment about ingredients necessary to create the self- 
organized critical state. This fact is reflected in the 
doubts about whether locally dissipative systems really 
present SOC or have only a very strong divergence of the 
mean avalanche size s when approaching the conserva- 
tive limit. The recent results by Chabanol and Hakin 0, 
Brock and Grassberger |3) and Kinouchi et al. Q stat- 
ing that the random-neighbor OFC model is not criti- 
cal in the dissipative regime and contradicting previous 
claims [g| , is a clear example of the difficulty of making 
such distinction solely on the basis of simulations. It is 
also worth remembering that the prototypical sandpilc 
(BTW) model is not critical in the presence of local dis- 
sipation 

The distinction between conservative/dissipative local 
dynamics, however, is not what is relevant for predict- 
ing critical behavior. The decisive ingredient seems to 
be the value of the coupling parameter J (or the nature 
of the distribution p{J) in non- homogeneous systems). 
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For example, the Feder and Feder model with k neigh- 
bors is non-conservative but is critical when the coupling 
constant is equal to J c = 1/k ||,|8| . 

In this paper, a model is proposed which is similar to, 
but simpler than, the random-neighbor stick-slip models 
studied in [0|j]. For this model, the stationary distribu- 
tion of states Poo(E) and the mean avalanche size s, as 
functions of the coupling parameter J, have simple an- 
alytical forms (in the limit of infinite system size) . The 
analysis in terms of branching processes is transparent 
and gives a clear mechanism for the emergence of very 
large but finite s in a non-negligible region of the param- 
eter space. In another words, although true criticality 
occurs only at a special point J c , there exist a large re- 
gion where power laws over several decades appear. In 
this region the behavior of the system can be considered 
almost critical. 

This occurs because the original parameter which con- 
trols the critical behavior (the branching rate a in a 
branching process) is now, in SOC models, a slow dy- 
namical variable cr t (J) that depends on the coupling pa- 
rameter J. In our model, the stationary value cr^,/) 
shows a plateau near the critical value a c = 1, thus en- 
larging the region in J space where the system displays a 
critical behavior. We will say that the system is critical 
for J = J c when a = 1 and is 'quasi-critical 'or 'almost 
critical ' for values of J where a ~ 1. This fact may be 
relevant as an explanation for the ubiquity of approxi- 
mate scale invariance in nature Q . 

The remainder of the paper is organized as follows: In 
Sec. II, the model is introduced and the main results 
obtained. The issue of robustness in SOC models is dis- 
cussed in Sec. III. Sec. IV contains concluding remarks 
and suggestions for future work. 

II. EXTREMAL FEDER AND FEDER MODEL 
(EFF MODEL) 

A. The model 

The EFF model is a random-neighbor version of the 
Feder and Feder model using an extremal dynamics 
similar to the Bak-Sneppcn model [^0[ . The extremal 
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dynamics, that in this case substitutes (and plays the 
same role that) the slow driving of the original Feder 
and Feder model, is here an essential ingredient for the 
observation of self-organized criticality. 

All sites j = 1, . . . , N have a continuous state variable 
Ej G 1Z. At each time step the site with maximal value 
'fires', resetting its value to zero plus a noise term rj. 
Then, k random 'neighbors' (rn) of the firing site have 
their states incremented by a constant J plus a noise 
term. The choice of neighbors is done at the firing in- 
stant: the randomness is annealed. So, denoting the ex- 
tremal value at instant t as E* = maxjEj}, the update 
rules are: 



»?(*). 



i) 

E rn {t + 1) = E rn (t) + J + r] rn (t), 



(1) 



with rj and rj rn being random variables uniformly dis- 
tributed in the interval [0, e] (the range of e will be dis- 
cussed later). Note that each random neighbor receives 
a different quantity r) rn . 

Consider the instantaneous density of states p t (E) . It 
is clear that for any E outside the intervals I n = [{n — 
1) J, (n — 1) J + ne],n = 1,2,..., this density decays to 
zero for long times. These intervals effectively discretize 
the phase space, so it is useful to define the following 
quantities, 



P, 



(n— 1) J+ne 



(n-l)J 



p(E) dE, 



(2) 



with n = 1, 2, . . . , n maxi and e < J/n max so that the in- 
tervals do not overlap (the integer n max will be obtained 
later). The process can be thought of as a transference 
of sites between the intervals /„ . At each time step, one 
site is transferred to the interval I\ and, with probability 
kPi, one site is removed from this interval. The aver- 
age flux to the intervals /„ with n > 1 corresponds to 
the probability kP n _\ that a neighbor is chosen in the 
previous I n -i interval minus the probability kP n that a 
neighbor is chosen in the interval /„. The average num- 
ber of sites in each interval is N n (t) — NP n (t). For long 
times, that is, when the density of states outside the /„ 
intervals goes to zero, one can write 



P 1 (t+l) = P 1 (t) + ±[l-kP 1 (t)}, 

P n {t + 1) = P n (t) + i [kP n _i(t) - kP n {t)] 



(3) 



Here, each time step is equal to the update of the maxi- 
mal site and k random neighbors. 

The condition for steady states, P n (t + 1) = P n (t) = 
P*, gives 



P* = l/k, 
p* _ p* 



(4) 



that is, P* = l/k for all n. But since p(E) is normalized, 
only n max intervals with P n of 0(1) can exist. That is, 



(5) 



giving that n max = fc.This means that Poo{E) is com- 
posed of k bumps (n = 1, . . . ,n rnax = k) and the pre- 
vious condition for producing non-overlapping intervals 
I n reads e < J/k. There is also a bump of 0(\ogN/N) 
(by analogy with the results from |ll|) situated at the 
interval Ik+i = [kJ, kJ + (k + l)e]. The other intervals 
n > k + 1 have P n of yet smaller order (see Fig. 1). 



B. Avalanches 

An avalanche will be defined as the number of firing 
sites until an extremal site value falls bellow the thresh- 
old E t h = 1 1 1 3 1 . Note that the first site of an avalanche 
(the 'seed') always has E < 1 but it counts as a fir- 



ing site. So, if a seed produces no supra-thrcshold sites 
('descendants'), this counts as an avalanche of size one. 
This definition of avalanches agrees with that used in the 
studies of relaxation oscillator models. 

In these random neighbor models, an avalanche can 
be identified as a branching process where an active site 
produces k new sites, each one having a probability p of 
being active (a 'branch') and a probability 1 — p of being 
inactive (a 'leaf'). The branching rate a = kp measures 
the probability that a firing site produces another firing 
site. 

A known result for a process with a constant branching 
rate a is the distribution of avalanche sizes 01 , 



(JN /CS-(S-I) 

k) 



(G) 



which, for large s and small S = 1 — a has the form 
1 



P(s) 



= S - 3 / 2 exp(-s/sA (7) 
V2t(1 - l/k) ^ ' * h 

2(*-l)„ _x-2 



k 



■(l-o-y 



(8) 



We will see that Eq. (jg) can be applied to the EFF model 
with the stationary value <Joo{J)- 

Now, consider an avalanche which has terminated after 
s sites have fired. This avalanche is composed of one 
seed and s — 1 descendants. But the average number of 
descendants produced by s firing sites is as. Thus, on 
average, the relation 



s — 1 = a s, 
must hold, which leads to 

1 



1 - cr 



(9) 



(10) 



Of course, this result can be obtained directly from 
Eq. (Q) after some work. Note that (Too = <r{t — > oo) 
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refers to the stationary value of the branching rate: dur- 
ing the transient, at changes with the avalanche time t. 
Although questioned by some authors ||, we retain the 
name self-organization for this evolution of a t toward Coo 
mainly as a label to distinguish these systems from stan- 
dard branching processes where a is fixed a 'priori. 



C. The J = l/k case 

In the case J = 1/fc, the calculation of a x is trivial. 
The fc-th bump, (n = n max ) which starts at (k — 1)J, 
must lie bellow the threshold E t h = 1 (if not, the system 
is supercritical). Then, e must satisfy the condition (fc — 
l)/k + fee < 1, that is, 



e < 1/fc 2 



(11) 



For the standard k = 4 neighbor case this reads e < 
0.0625. This condition also implies that neighbors per- 
taining to the other bumps do not contribute to a^ , that 
is, cannot fire when receiving a maximal contribution 
J + e. Now, since all the neighbors pertaining to the 
fc-th bump receive at least the quantity J = 1/fc, they 
are always transformed into active sites. Thus, the aver- 
age number of descendants of a firing site is 



kPt = k x 



1 



1. 



(12) 



which corresponds to a critical branching process. It is 
known that in this case the system presents an infinite s 
(see Eq. (|l0|)) and, for large s, a pure power law 



P{s) = 



1 



-3/2 



(13) 



v/2tt(1 - 1/fc) 
for the distribution of avalanche sizes j|. 

D. Results for general J 

For the case J < 1/fc, in order to obtain an expression 
for Coo (J), the knowledge of the distribution of states 
Poo(E) is required. But it is clear that if kJ = 1 — 6 
then inevitably er^ < 1 (even for very small S > 0), 
since some sites pertaining to the fc-th bump may not 
receive a sufficient contribution to make them active (see 
Eq. (??) below). Thus, any value J < J c = 1/fc is sub- 
critical. This is a common feature of many models with 

socna. 

In our model, the calculation of p oa (E) is very simple. 
In the stationary state, a site pertaining to the n-th bump 
has energy E = (n — 1) J + z n , where z n is the sum of 
n random variables uniformly distributed in the interval 
[0, e]. The distribution p(z n ) may be calculated from 

p( Zl ) -e- 1 9(z 1 )e(e-z 1 ), 

/>OC 

P(z n +l) = / dz n dzip(z n )p(zi)S(z n + Z\ - Z n+1 ) . 



For the k = 4 case, 

p(z 2 ) = e~ 2 z 2 Q(z 2 ) 6(e - z 2 ) 

+ (2e-z 2 )e(z 2 -e)6(2e-z 2 ), 

p(z 3 ) = r 3 |e(z 3 )e( f - Z 3) 

3e 2 



(14) 



-z 2 + 3ez 3 - — 6(z 3 - e) 9(2e - z 3 ) 



z 3 ~ 9e 



y - 3ez 3 + — ) 6(z 3 - 2e) 6(3e - z 3 ), 



p(z A ) = e - 4 -^-e(z 4 )e(e - z 4 ) 



2e 3 
3 

2e 3 



^ + 2ez 2 - 2e 2 z 4 + ^ ) 6(z 4 - e)9(2e - z 4 ) 



— + 2e.x 2 - 2e 2 a; + ^— ) 6(z 4 - 2e)9(3e - z 4 ) 



+ — 9(z 4 - 3e)9(4e - z 4 ), 
6 



(15) 



with the shorthand x = (4e — z 4 ). The distribution 
Poo^E) has k bumps. Each bump (labeled by n) starts at 
E n = (n — 1) J, being proportional to p(z n ) (the constant 
of proportionality is just 1/fc). In Fig. 1, the distribution 
Poo(E) is compared with simulation results for a system 
with 10 sites, J = 0.235, e = 0.05 and a sufficient num- 
ber of avalanches. 

For such large systems, we must be careful about using 
reliable random neighbor generators. In order to speed 
up the search for the extremal site, we used the binary 
rooted tree algorithm described by Grassberger [jl2j . For 
example, if the system has 2 m sites, a binary tree with 
m+1 levels is created such that, in each node at level /, it 
is stored the largest value of E of the two branch nodes 
of the I + 1-th level. So, the 0-th (root) level contains 
the value of the extremal site. Ascending the tree, we 
locate the position of this site in the upper level. After 
the extremal site firing, the tree must be updated. The 
same occurs when the random neighbors are updated. 
These operations have a complexity 0(\ogJV) instead of 
the 0(Af) complexity of the naive search mechanism. 

The stationary branching rate CTqc is calculated as fol- 
lows. All the sites that can be activated pertain to the fc- 
th bump. When hit, sites with E > 1 — J are always acti- 
vated. In terms of the re-scaled variable 2^ = E—(k—l)J, 
this condition refers to sites with Zk > 5 = 1 — k J. They 
contribute to the branching rate with the quantity a', 



p(E) dE 



i-j 



s+j 



p(z)dz, 



(16) 



where z = Zk- 

Sites with E < 1 — J — e cannot be activated and do 
not contribute to a. Sites with 1 — J — e < E < 1 — J 
can be activated if they receive a quantity J + r/ > 1 — E, 
that is, t] > 5 — z. This occurs with probability P(n > 
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6 — z) = 1 — (6 — z)/e. Thus, these sites contribute to the 
branching rate with the quantity 

/■lkJ 

a" = k / P(E) P(r) > 1 - E - J) dE, 

Jl-kJ-e 
r<5 



p(z) ( 1 ) </: . 

S-e V e 



The total branching rate is then 

r S-e 



(17) 



Co 



a' + a" 



1 



p(z) dz 



S 1 f'^ 

— / p(z) dzH — / zp(z)dz, 

e e J s-e 



(18) 



where we used the fact that J^ +J p(z)dz — 1. Since p(z) 
has a simple piece- wise polynomial form (see Eq. (|l4|)) 
the 

calculation of s is straightforward and the result is pre- 
sented in Fig. 2 along with simulation results for the 
k — 4, e = 0.05, for systems with up to N = 2 18 = 
262 144 sites. In Fig. 3, we plot simulation results for the 
P(s) distribution which agree very well with Eq. (0) if 
a = (?ao{J) is used in that expression. Strong finite size 
effects, however, are present when J > 0.235. 

For 5 < e, that is, J c — J < e/k, the form assumed by 
CToo is particularly simple, since p(z) = Ce~ k z k ~ 1 in that 
interval (C is a numerical constant). Then, 



Coo = 1 — 



= 1 - 



C 

e k+1 jo 
C 

k(k + l) 



z k ~\b-z)dz 



fc+1 



(19) 



the avalanche cutoff lenght 

Since 6 = 1 — kJ = k(J c — J), we obtain, from 
Eqs.( ||) and (0), the avalanche cutoff size and the aver- 
age avalanche size 

_ 2(fc-l)(fc + l)V(^) 
s e- c^fc+i (Jc-J) , (20) 

(fc+l)e^ /2 

with the critical exponent 

u = 2(k + l). (21) 

For example, with k = 4 (which means C = 1/6, see 
Eq. @) and e = 0.05, the mean avalanche size is s = 120 
already for J = 0.2375. Curiously, this behavior is similar 
to the s oc (J c — J)~ k divergence found in the standard 
random neighbor FF model |3J. 



E. The EFF model with noiseless couplings 

It is instructive to compare the above behavior with 
that of a simpler EFF model pj where the firing rule is 
the same, E*(t+1) = r\ S [0, e], but the coupling between 
sites is noiseless, E rn (t + 1) = E rn {t) + J. Thus, Poc{E) 
assumes the form of k rectangular bumps with p(z n ) — 
e~ 1 0(z n )0(e — z n ). In this noiseless EFF model, the 
branching rate, the cutoff size and the average avalanche 
size are 







for 6 > e 



1-6/e for < 6 < e ' 



^ = 2(fc-l)p(J c - J)- 2 



k 



(Jc - jy 1 



(22) 



In contrast with the noisy model, large avalanches only 
occur when J is very close to J c (see Fig. 2). Thus, the 
EFF model with noiseless couplings does not present an 
enlargement of the region where the system displays a 
critical behavior as observed in the noisy EFF model. 



III. ON SOC DEFINITIONS 

The idea of self-organized criticality present in the lit- 
erature embodies two distinct properties. The term crit- 
ical refers to the existence of power laws and to the ab- 
sence of a characteristic scale in the response of the sys- 
tem to the driving mechanism of the dynamics; the term 
self-organized refers to the fact that there exist a param- 
eter (at), which controls the avalanching process, whose 
value is not fixed a priori like, for example, in standard 
percolation and branching processes. This parameter 
evolves in time, during a transient phase, toward a sta- 
tionary value (Too. Indeed, this time dependence should 
be written as at = a(pt(E)), that is, a t is a functional 
of the distribution of states pt(E), that, in turn, evolves 
toward a statistically stationary distribution p oa (E). So, 
CToo = a(poo(E)). If CToo = a c = 1, the system is critical. 

The evolution of p t (E) toward the steady-state Poo (E) 
is akin to the transient relaxation in equilibrium systems: 
any initial condition leads to the same stationary state, 
thus to the same value of CToo- However, this robustness 
to initial conditions and external perturbations on p(E) 
('dynamical stability') should not be mistaken as param- 
eter robustness ('structural stability'). This is a distinct 
characteristic claimed to be present on some SOC mod- 
els (see for instance ||,|l^,|l^,|l8| ) . For a system to have 
'structural stable criticality', there would be a finite pa- 
rameter range for which, after the transient, the system 
is critical. In this CH.SC, (J qq (J) = a c for J belonging to 
some interval [ J c , 1/k] . That kind of behavior will be also 
called by us generic SOC . 

'Structural stability' is a relative concept which de- 
pends on the parameter space physically available for the 
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system. For example, it is well known that the sandpile 
model is not critical in the presence of dissipation. The 
sandpile dissipation parameter corresponds to the quan- 
tity 6 = 1 - kJ in our model |,||. The standard BTW 
model is by definition 'tuned' into a critical state through 
the 'imposition' of a conservation law. Although it could 
be argued that dissipation is not a natural feature of 
sandpiles, since sand does not disappear, the appearence 
of SOC in nature would sound much more natural if crit- 
icality could be observed over a region of the parameter 
space, not only in a special point. 

Generic self-organized criticality is depicted in 
curve (a) of figure 4. In this case, there is a finite range of 
J values for which assumes the critical value <r c = 1 . 
In this figure, curves (d) and (e) represent the behav- 
ior observed in the BTW model and also in the noise- 
less EFF model examined above, for which the system 
is critical only for a special value of the parameter J. 
However, there is a third possibility. Curves (b) and (c) 
represent the behavior of a (J) given by Eq. ([18|) for the 
EFF model with noisy couplings: although the system 
is critical only at J = J c , the system is 'almost critical' 
over a large parameter region. This behavior has also 
been observed in the standard random neighbor versions 
of FF and OFC models || . The importance of this char- 
acterization is that several models in the SOC literature, 
previously seen as having true generic criticality, are now 
recognized as having only an almost critical behavior as 
discussed above. 

A model which apparently presents generic SOC be- 
havior in coupling space is the two-dimensional OFC 
model | il||l7| , |l9f . Also the standard Feder and Feder 
model [|| is claimed to be critical for J < J c p8| , |l9| | . 
Looking at the behavior of the models studied so far, 
we make the following conjecture: a necessary condition 
for a lattice model to present a generic SOC behavior is 
that its corresponding random neighbor version already 
presents an enlarged critical region in the sense discussed 
above. This could be tested by comparing the 2D ver- 
sions of the EFF and noiseless EFF models studied above. 

In conclusion, we found that some systems that display 
SOC, although being critical only for a single value for 
J, are almost critical region in a large region of the pa- 
rameter space. This almost critical behavior is difficult 
to be distinguished, in practice, from true generic SOC 
behavior: both in numerical simulations (huge lattices 
would have to be used) and in Nature (due to limitations 
in the data) power laws can only be measured over some 
scale decades M. So, in order to explain the ubiquity of 
scale invariance in Nature, having a true generic SOC or 
only presenting an enlarged region where the system is 
almost critical are, as far one can measure, identical. 



IV. CONCLUSIONS 

A class of extremal stick-slip models has been intro- 
duced and studied in the N — > oo limit. We showed that 
noise in the couplings of the EFF model changes the ex- 
ponent that controls the amplitude of the critical region 
from v = 2 to v = 2(k + 1). This enlargement of the 
region where the system displays a critical behavior is 
similar to that found in the standard random neighbor 
OFC and FF models [||-[|]. Like in other models, the 
true critical state occurs only for one point in parameter 
space P,p|,p|j7p^| , but in practice that fact can hardly be 
noticed, and the model displays the typical features of 
generic SOC. 

In future work we hope to determine the minimal in- 
gredients for producing the enlargement of the critical 
region in the models examined in the SOC literature. We 
will also present results for the two-dimensional case and 
compare with the standard OFC and FF models. The 
simple mechanism devised in this work suggests that, if 
true generic criticality is not easy to obtain in the space 
of possible models, this quasi-critical behavior certainly 
is. Thus, for explaining the robustness of approximate 
scale invariance in Nature, this mechanism seems to be 
more "generic" than generic criticality. 
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FIGURE CAPTIONS 

Figure 1: Distribution of states Poo (E) for k = 4, J = 
0.235 and e = 0.05: theoretical (solid) and simulation 
(circles) with N = 10 4 sites. 

Figure 2: Mean avalanche size s as a function of pa- 
rameter J. Theoretical (solid) and simulations with N 
up to 2 18 = 262 144 sites for noisy EFF model (circles) 
with e = 0.05 and noiseless EFF model(triangles) with 
e = 0.2. These e values are chosen such that the last 
interval (I 4) has the same length in both models. 

Figure 3: Simulation results (N = 2 13 = 8192 
sites, k — 4, e = 0.05) for the distribution P(s) with 
J = 0.21,0.22,0.23,0.235 (from left to right), compared 
with theoretical curves (solid). 

Figure 4: a) Generic self-organized criticality: the 
value of parameter Coo is critical on a finite range of the 
system parameter J; b) e = 0.0625 and c) e = 0.05, 
enlargement of the critical region (EFF model with noisy 
couplings, k = 4): is almost constant near J c ; d) e = 
0.25 = 4 x 0.0625 and e) e = 0.2 = 4 x 0.05, standard 
critical behavior (EFF model with noiseless couplings): 
the coupling parameter J must be very close to 0.25 for 
obtaining er^ w a c due to the linear behavior of (Too (J). 
Note that, in the noiseless couplings case, e refers to the 
amplitude of the noise received by the extremal site after 
discharge. 
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Fig. 1 - Kinouchi and Prado 




Fig. 2 - Kinouchi and Prado 
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Fig. 4 - Kinouchi and Prado 
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